READ ME

This R code is used to take raw respirometry slopes measured on Galaxias maculatus at Deakin’s Marine Research Centre in Queenscliff, Victoria, Australia and calculate their standard, routine, and maximum metabolic rate. These metabolic measurements are made repeatedly on known individuals held under one of two temperature treatments (18°C or 23°C), and paired with repeated mass and length measurements over five month of data collection (May - September, 2023). After reading in all data, a multivariate mixed effects modeling approach is undertaken in order to assess covariance between growth and metabolism traits among- and within-individuals.

LOAD REQUIRED PACKAGES

# installs and loads packages, need to install "pacman" first

pacman::p_load("brms",
               "broom",
               "broom.mixed",
               "car",
               "coda",
               "dplyr",
               "ggeffects",
               "ggiraphExtra",
               "ggplot2",
               "ggpubr",
               "hms",
               "lme4",
               "lmerTest",
               "lubridate",
               "nlme",
               "posterior",
               "readxl",
               "tidyverse"
)

Working directories

wd <- getwd() # getwd tells us what the current wd is, we are using this to drop it in a variable called wd

raw_data_wd <- paste0(wd, "./00_raw_data")  # creates a variable with the name of the wd where raw data are stored

    MRcalcs_wd <- paste0(raw_data_wd, "./MRcalcs")
    
    MR_slopes_wd <- paste0(raw_data_wd, "./MR_Slopes")
    
    MMR_wd <- paste0(raw_data_wd, "./MMR")

metadata_wd <- paste0(wd, "./01_metadata")  # creates a variable with the name of the wd where metadata are stored

figures_wd <- paste0(wd, "./03_figures")  # creates a variable with the name of the wd where figures are to be stored

Read in data, calculate metabolic rates, merge into master data table

#calcSMR function adapted from Chabot, 2016 (https://doi.org/10.1111/jfb.12845)
#note that we have not used the mlnd function as the mean of the lowest 10th percentile was better suited to our data

calcSMR <- function(Y) {
  u <- sort(Y)
  tenpc <- round(0.1 * length(u))
  SD10pc <- sd(u[1:tenpc])
  low10pc = mean(u[(which((u > (mean(u[1:tenpc])-SD10pc)))):(tenpc+which((u > (mean(u[1:tenpc])-SD10pc -u[1]))))])
  return(list(low10pc = low10pc))
}

# Define a function to extract date strings from file names
extract_date <- function(files) {
  gsub(".*?([0-9]{8}).*", "\\1", basename(files))
}

# Specify the directories where your files are located
csv_folder <- MR_slopes_wd
xlsx_folder <- MMR_wd

# Get the list of files in each folder
csv_files <- list.files(csv_folder, pattern = "_MR_slopes.csv", full.names = TRUE)
xlsx_files <- list.files(xlsx_folder, pattern = "_MMR.xlsx", full.names = TRUE)


# Extract date strings from file names
csv_dates <- extract_date(csv_files)
xlsx_dates <- extract_date(xlsx_files)


matching_files <- Map(function(date) {
  csv_file <- csv_files[grep(date, csv_dates)]
  xlsx_file <- xlsx_files[grep(date, xlsx_dates)]
  return(list(csv = csv_file, xlsx = xlsx_file))
}, unique(c(csv_dates, xlsx_dates)))


# Loop through matching files and perform your data processing
for (i in seq_along(matching_files)) {
  date <- unique(csv_dates)[i]
  current_files <- matching_files[[i]]  # Renamed 'files' to 'current_files'
  
  # Process CSV data
  tb_respirometry <- read_csv(current_files$csv) %>%
    rename(
      chamber_ch = Chamber.No,
      ID_fish    = Ind,
      mass       = Mass,
      length     = Length,
      volume_ch  = Ch.Volume,
      DOunit     = DO.unit,
      dateTime   = Date.Time,
      Temp_class = Temp.class,
      slope_wBR  = Slope.with.BR,
      BRSlope    = BR.Slope
    )  %>%
    mutate(
      dateTime       = as.POSIXct(ymd_hms(dateTime)),
      Time           = as_hms(ymd_hms(dateTime)),
      Date           = as.Date(dateTime, format = "%Y/%m/%d"),
    )
  
  
  tb_respirometry <- tb_respirometry %>%
    mutate(
      volume_net = volume_ch - mass,
      MR_wBR     = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
      BR         = BRSlope*(volume_ch/1000)*60*60,
      MR         = MR_wBR + BR,
      BR_perc    = (BR/MR_wBR)*100
    ) 
  
  #Calculate RMR and variance in MR for each fish, create a new table for individual RMR calcs 
  
  tb_rmr <- tb_respirometry %>%
    group_by(chamber_ch) %>%
    arrange(chamber_ch, dateTime) %>%
    slice(3:(n() - 1)) %>%
    ungroup()  %>%
    group_by(
      ID_fish, mass, length, volume_net, chamber_ch, Temp_class
    ) %>% 
    arrange(ID_fish) %>%
    drop_na() %>%
    summarise(
      RMR        = mean(MR),
      RMR_var    = var(MR),
      RMR_perc   = (RMR_var/RMR)*100,
      time_start = dateTime %>% min(),
      time_end   = dateTime %>% max()
    ) 
  
  #Calculate SMR for each fish, create a new table for individual SMR calcs 
  
  tb_smr <-
    tb_respirometry %>%
    group_by(
      ID_fish, mass,length, volume_net, chamber_ch
    ) %>% 
    arrange(ID_fish) %>%
    drop_na() %>%
    summarise(
      SMR        = calcSMR(MR)$low10pc %>% unname(),
      time_start = dateTime %>% min(),
      time_end   = dateTime %>% max()
    ) 
  
    # Process XLSX data
    tb_mmr <- read_excel(current_files$xlsx) %>%
      rename(
        chamber_ch = Chamber.No,
        ID_fish    = Ind,
        mass       = Mass,
        volume_ch  = Ch.Volume,
        DOunit     = DO.unit,
        dateTime   = Date.Time,
        Temp_class = Temp.class,
        slope_wBR  = Slope.with.BR,
        BRSlope    = BR.Slope
      )  %>%
      mutate(
        dateTime       = as.POSIXct(ymd_hms(dateTime)),
        Time           = as_hms(ymd_hms(dateTime)),
        Date           = as.Date(dateTime, format = "%Y/%m/%d"),
      ) %>%
      drop_na()
    
    tb_mmr <- tb_mmr %>%
      mutate(
        volume_net = volume_ch - mass,
        MMR_wBR    = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
        BR         = abs(BRSlope)*(volume_ch/1000)*60*60,
        MMR        = MMR_wBR - BR,
        BR_perc    = (BR/MMR_wBR)*100
      ) %>%
      arrange(ID_fish)
      
    # Combine the processed data
    # Combine the processed data, selecting only specific columns
    tb_MR_master <- tb_rmr %>%
      left_join(select(tb_smr, ID_fish, SMR), by = "ID_fish") %>%
      left_join(select(tb_mmr, ID_fish, MMR), by = "ID_fish") %>%
      select(-matches(".x$")) %>%
      rename_with(~gsub("\\.y$", "", .), matches(".y$"))
    
  # Export data
  setwd(MRcalcs_wd)
  write.csv(tb_MR_master, file = paste0(date, "_MRcalcs.csv"), col.names = NA, row.names = FALSE)
}

# read all "_MRcalcs.csv" files, bind into one table, and add columns with logged values 
file.list <- list.files(MRcalcs_wd)

setwd(MRcalcs_wd)
calcs_all <- lapply(file.list, read_csv)

tb_MRcalcs <- bind_rows(calcs_all, .id="Resp_Day") %>%
  mutate(
    log_SMR_low10pc = log10(SMR),
    log_mass = log10(mass),
    log_RMR = log10(RMR),
    log_MMR = log10(MMR),
    Resp_Day = as.numeric(Resp_Day)
  )

# assign Month category to all timepoints 
tb_MRcalcs$Month <- cut(
  tb_MRcalcs$Resp_Day,
  breaks = c(0, 13, seq(23, max(tb_MRcalcs$Resp_Day) + 10, by = 10)),
  right = FALSE,
  labels = c("April", "May", "June", "July", "August", "September")
)

setwd(raw_data_wd)
biometrics <- read.csv("FinalBiometrics.csv")

ds <- left_join(tb_MRcalcs, biometrics, by = c("ID_fish")) #a table with all resp data and metadata collected in dissection

ds1 <- ds[which(ds$Month != "April"),] # remove April data from analysis due to temperature inconsistency

Multivariate mixed-effects models to test the relationships between growth and metabolism in Galaxias maculatus.

ds1 <- ds1 %>%
  # Filter to keep only ID_fish where a mass exists for Month == "May"
  filter(ID_fish %in% ID_fish[Month == "May" & !is.na(mass)])


month_order <- c("May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)
ds1$Month <- factor(ds1$Month, levels = month_order) # assign month_order to factor variable "Month" in ds1

min_mass = min(ds1$mass)
min_log_mass = min(ds1$log_mass)

min_finalmass = min(((ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
min_log_finalmass = min((log10(ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))


# format data columns and generate centred and log-scaled variables
ds1 <- ds1 %>%
  mutate(
    mass_centre = mass - min_mass,                                  # centred mass  
    log_mass_centre = log_mass - log10(min_mass),                   # log-scaled centred mass  
    Temp_class = as.factor(Temp_class),                             # set temperature treatments as factors
    SMR1 = log10(SMR),                                              # log-scaled SMR
    RMR1 = log10(RMR),                                              # log-scaled RMR
    MMR1 = log10(MMR),                                              # log-scaled MMR
    mass1 = log_mass,                                               # log-scaled mass
    Month_int = case_when(Month == "May" ~ 1,
                          Month == "June" ~ 2,
                          Month == "July" ~ 3,
                          Month == "August" ~ 4,
                          Month == "September" ~ 5),
    Month_centred = Month_int - 1,                                  # left-centred mass so May is the intercept
    Month_centred_factor = as.factor(Month_centred),                # set centred month as factors
    Rack = as.integer(Rack),                                        # set rack ID to integer (for now)
    Bin = as.integer(Bin),                                          # set bin ID to integer (for now)
    finalmass_centre = Final_Mass_g - min_finalmass,                # centre last mass measurement taken
    log_finalmass_centre = log10(Final_Mass_g) - min_log_finalmass, # log-scaled and centred final mass
    log_fatmass = log10(Fat_mass_g)                                 # log-scaled fat mass
  )

ds1 <- tidyr::unite(data = ds1, Rack, Bin, col = "tankID", sep = "_") # merge rack and bin data into one identifier, call it tankID

#read in density data table, call it BinCounts
setwd(raw_data_wd)
BinCounts <- read_excel("BinCounts.xlsx", sheet = "data")

#merge BinCounts and ds1
ds1 <- left_join(ds1, BinCounts %>% dplyr::select("Month", "tankID", "Population"), by = c("Month", "tankID"))

#reformat Sex variable so that both unknown and immature fish are "NA" 
ds1$Sex <- ifelse(ds1$Sex == "I", NA, ds1$Sex)

# rename "population" as density, set tankID as factor
ds1 <- ds1 %>%
  rename(
    Density = Population
  ) %>%
  mutate(
    tankID = as.factor(tankID)
  )

# check structure of ds1
str(ds1)
## tibble [624 × 41] (S3: tbl_df/tbl/data.frame)
##  $ Resp_Day            : num [1:624] 13 13 13 13 13 13 13 13 13 13 ...
##  $ ID_fish             : num [1:624] 81 82 83 84 85 86 87 88 89 90 ...
##  $ chamber_ch          : chr [1:624] "A1" "B3" "B4" "A2" ...
##  $ Temp_class          : Factor w/ 2 levels "18","23": 2 2 2 2 2 2 2 2 2 2 ...
##  $ RMR                 : num [1:624] 0.952 0.412 0.712 1.016 0.876 ...
##  $ RMR_var             : num [1:624] 0.0149 0.0135 0.0429 0.0696 0.0454 ...
##  $ RMR_perc            : num [1:624] 1.57 3.28 6.02 6.85 5.18 ...
##  $ time_start          : POSIXct[1:624], format: "2023-05-22 17:53:00" "2023-05-22 17:53:00" ...
##  $ time_end            : POSIXct[1:624], format: "2023-05-23 08:03:00" "2023-05-23 08:03:00" ...
##  $ mass                : num [1:624] 2.65 0.81 1.78 2.73 2.29 0.98 2.04 2.43 2.31 1.65 ...
##  $ length              : num [1:624] 77 62 73 78 80 64 80 77 78 68 ...
##  $ volume_net          : num [1:624] 707 709 708 707 708 ...
##  $ SMR                 : num [1:624] 0.845 0.284 0.53 0.804 0.707 ...
##  $ MMR                 : num [1:624] 1.43 1.11 1.44 1.64 1.69 ...
##  $ log_SMR_low10pc     : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
##  $ log_mass            : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
##  $ log_RMR             : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
##  $ log_MMR             : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
##  $ Month               : chr [1:624] "May" "May" "May" "May" ...
##  $ tankID              : Factor w/ 20 levels "1_10","1_11",..: 17 17 17 17 17 17 17 17 18 18 ...
##  $ TagID               : chr [1:624] "YL" "YR" "BrL" "BrR" ...
##  $ Final_Mass_g        : num [1:624] 4.62 3.65 2.64 6.33 5.19 ...
##  $ Final_Length_mm     : int [1:624] 94 86 84 98 100 95 90 89 91 98 ...
##  $ Gonad_mass_g        : num [1:624] 0.035 0.022 0.018 0.039 0.052 0.028 0.362 0.033 0.015 0.093 ...
##  $ Sex                 : chr [1:624] "F" "M" "F" "M" ...
##  $ Fat_mass_g          : num [1:624] 0.315 0.123 0.062 0.314 0.214 0.091 0.076 0.044 0.15 0.209 ...
##  $ Ventricle_mass_g    : chr [1:624] "0.017" "" "" "" ...
##  $ Otolith_age_y       : num [1:624] 1 NA NA NA NA 1 1 NA NA 1 ...
##  $ mass_centre         : num [1:624] 2.07 0.23 1.2 2.15 1.71 0.4 1.46 1.85 1.73 1.07 ...
##  $ log_mass_centre     : num [1:624] 0.66 0.145 0.487 0.673 0.596 ...
##  $ SMR1                : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
##  $ RMR1                : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
##  $ MMR1                : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
##  $ mass1               : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
##  $ Month_int           : num [1:624] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Month_centred       : num [1:624] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Month_centred_factor: Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ finalmass_centre    : num [1:624] 3.3 2.33 1.32 5.01 3.87 ...
##  $ log_finalmass_centre: num [1:624] 0.544 0.442 0.301 0.681 0.594 ...
##  $ log_fatmass         : num [1:624] -0.502 -0.91 -1.208 -0.503 -0.67 ...
##  $ Density             : num [1:624] 8 8 8 8 8 8 8 8 5 5 ...
head(ds1)
## # A tibble: 6 × 41
##   Resp_Day ID_fish chamber_ch Temp_class   RMR RMR_var RMR_perc
##      <dbl>   <dbl> <chr>      <fct>      <dbl>   <dbl>    <dbl>
## 1       13      81 A1         23         0.952  0.0149     1.57
## 2       13      82 B3         23         0.412  0.0135     3.28
## 3       13      83 B4         23         0.712  0.0429     6.02
## 4       13      84 A2         23         1.02   0.0696     6.85
## 5       13      85 A4         23         0.876  0.0454     5.18
## 6       13      86 A3         23         0.751  0.0504     6.71
## # ℹ 34 more variables: time_start <dttm>, time_end <dttm>, mass <dbl>,
## #   length <dbl>, volume_net <dbl>, SMR <dbl>, MMR <dbl>,
## #   log_SMR_low10pc <dbl>, log_mass <dbl>, log_RMR <dbl>, log_MMR <dbl>,
## #   Month <chr>, tankID <fct>, TagID <chr>, Final_Mass_g <dbl>,
## #   Final_Length_mm <int>, Gonad_mass_g <dbl>, Sex <chr>, Fat_mass_g <dbl>,
## #   Ventricle_mass_g <chr>, Otolith_age_y <dbl>, mass_centre <dbl>,
## #   log_mass_centre <dbl>, SMR1 <dbl>, RMR1 <dbl>, MMR1 <dbl>, mass1 <dbl>, …

GROWTH ANALYSES

Visualize growth data

hist(ds1$mass) # plot all mass values to see mass range and most frequent values

ds1$Month <- factor(ds1$Month, levels = month_order)              # reassign month_order to factor variable "Month" in ds1

# Plot mass over time, connect by ID_fish to show individual trends

ggplot(ds1, aes(x = Month, y = mass, group = ID_fish)) +          # plot raw data by ID
  geom_point(size = 2) +                                          # size of data points
  geom_line() +                                                   # connect points by individual ID
  facet_wrap(. ~ Temp_class) +                                    # produces one plot per temperature treatment
  theme_bw()

Log-linearize mass data

We can also visualise this by plotting log10 mass over log10_x (note that we cannot use log10_y as some masses were < 1 initially)

ggplot(ds1, aes(x = Month_int, y = log_mass, group = ID_fish)) +  # plot logged mass data by ID
  geom_point(size = 2) +                                          # size of data points
  geom_line() +                                                   # connect points by individual ID
  facet_wrap(. ~ Temp_class) +                                    # produces one plot per temperature treatment
  theme_bw() +
  scale_x_log10()                                                 # display x-axis on log scale

Univariate growth models for two temperature treatments

ds18 <- ds1[which(ds1$Temp_class == "18"),]                   # subset ds for only fish at 18C
ds23 <- ds1[which(ds1$Temp_class == "23"),]                   # subset ds for only fish at 23C

mass18_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 18C only
           data=ds18, REML = T)

summary(mass18_mod)                                 # generate model output
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -492
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8214 -0.4648  0.0139  0.4865  2.3450 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.040770 0.20191       
##           Month_int   0.002038 0.04515  -0.57
##  Residual             0.003528 0.05939       
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)     0.32789    0.02232  14.689
## MonthJune       0.15461    0.01147  13.478
## MonthJuly       0.22113    0.01498  14.761
## MonthAugust     0.15515    0.01989   7.799
## MonthSeptember  0.36370    0.02487  14.622
## 
## Correlation of Fixed Effects:
##             (Intr) MnthJn MnthJl MnthAg
## MonthJune   -0.364                     
## MonthJuly   -0.412  0.629              
## MonthAugust -0.410  0.600  0.770       
## MonthSptmbr -0.408  0.581  0.776  0.865
plot(mass18_mod)                                    # residuals plot

mass23_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 23C only
                 data=ds23, REML=T)

summary(mass23_mod)                                 # generate model output
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -515.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9289 -0.3941 -0.0018  0.5057  4.8187 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.040995 0.20247       
##           Month_int   0.001445 0.03802  -0.46
##  Residual             0.003138 0.05602       
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)     0.21996    0.02451   8.974
## MonthJune       0.12091    0.01123  10.771
## MonthJuly       0.21448    0.01402  15.297
## MonthAugust     0.22778    0.01773  12.848
## MonthSeptember  0.35157    0.02193  16.028
## 
## Correlation of Fixed Effects:
##             (Intr) MnthJn MnthJl MnthAg
## MonthJune   -0.298                     
## MonthJuly   -0.335  0.623              
## MonthAugust -0.341  0.611  0.773       
## MonthSptmbr -0.337  0.589  0.778  0.858
plot(mass23_mod)                                    # residuals plot

Combined growth model

# mass model including both temperature treatments, and their interaction with Month
mass_mod <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish), 
                 data=ds1, REML=T)

summary(mass_mod)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 +  
##     Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -789
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6831 -0.4303  0.0066  0.4509  4.5846 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.038839 0.19708       
##           Month_int   0.001745 0.04177  -0.55
##  Residual             0.003485 0.05903       
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                              0.175461   0.059281 337.150651   2.960
## MonthJune                                0.156983   0.011479 391.676251  13.676
## MonthJuly                                0.224031   0.014800 247.773510  15.137
## MonthAugust                              0.162564   0.020135 174.402058   8.074
## MonthSeptember                           0.373693   0.024476 123.702938  15.268
## as.factor(Temp_class)23                 -0.071983   0.036905 120.004217  -1.950
## SexF                                     0.123289   0.043645 115.959487   2.825
## SexM                                     0.051217   0.046603 115.166678   1.099
## Density                                  0.010421   0.005792 468.638813   1.799
## MonthJune:as.factor(Temp_class)23       -0.026407   0.018735 390.916591  -1.409
## MonthJuly:as.factor(Temp_class)23       -0.003931   0.023841 236.415811  -0.165
## MonthAugust:as.factor(Temp_class)23      0.065838   0.031139 146.333633   2.114
## MonthSeptember:as.factor(Temp_class)23  -0.015746   0.038370 107.376799  -0.410
##                                        Pr(>|t|)    
## (Intercept)                             0.00330 ** 
## MonthJune                               < 2e-16 ***
## MonthJuly                               < 2e-16 ***
## MonthAugust                            1.07e-13 ***
## MonthSeptember                          < 2e-16 ***
## as.factor(Temp_class)23                 0.05345 .  
## SexF                                    0.00557 ** 
## SexM                                    0.27405    
## Density                                 0.07265 .  
## MonthJune:as.factor(Temp_class)23       0.15948    
## MonthJuly:as.factor(Temp_class)23       0.86916    
## MonthAugust:as.factor(Temp_class)23     0.03618 *  
## MonthSeptember:as.factor(Temp_class)23  0.68236    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mass_mod)      # residuals plot

ranova(mass_mod)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + Month:as.factor(Temp_class)
##                                        npar logLik     AIC   LRT Df Pr(>Chisq)
## <none>                                   17 394.50 -755.01                    
## Month_int in (1 + Month_int | ID_fish)   15 331.56 -633.11 125.9  2  < 2.2e-16
##                                           
## <none>                                    
## Month_int in (1 + Month_int | ID_fish) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mass model including both temperature treatments, and their interaction with Month
mass_mod_noMonth <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1|ID_fish), 
                 data=ds1, REML=T)

summary(mass_mod_noMonth)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 |  
##     ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -663.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6703 -0.5089  0.0290  0.4797  3.9785 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.026454 0.16265 
##  Residual             0.007503 0.08662 
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                              0.193236   0.060932 357.632744   3.171
## MonthJune                                0.156704   0.015013 391.702422  10.438
## MonthJuly                                0.221347   0.015486 395.495181  14.293
## MonthAugust                              0.161419   0.017561 412.738321   9.192
## MonthSeptember                           0.373567   0.017561 412.738321  21.273
## as.factor(Temp_class)23                 -0.073410   0.036398 155.881117  -2.017
## SexF                                     0.121574   0.043722 118.677173   2.781
## SexM                                     0.050315   0.046686 116.817268   1.078
## Density                                  0.008195   0.006125 491.655639   1.338
## MonthJune:as.factor(Temp_class)23       -0.026128   0.024591 391.220473  -1.062
## MonthJuly:as.factor(Temp_class)23       -0.001748   0.024840 392.469865  -0.070
## MonthAugust:as.factor(Temp_class)23      0.065981   0.026055 399.871498   2.532
## MonthSeptember:as.factor(Temp_class)23  -0.015885   0.026056 398.802053  -0.610
##                                        Pr(>|t|)    
## (Intercept)                             0.00165 ** 
## MonthJune                               < 2e-16 ***
## MonthJuly                               < 2e-16 ***
## MonthAugust                             < 2e-16 ***
## MonthSeptember                          < 2e-16 ***
## as.factor(Temp_class)23                 0.04543 *  
## SexF                                    0.00631 ** 
## SexM                                    0.28337    
## Density                                 0.18154    
## MonthJune:as.factor(Temp_class)23       0.28868    
## MonthJuly:as.factor(Temp_class)23       0.94394    
## MonthAugust:as.factor(Temp_class)23     0.01171 *  
## MonthSeptember:as.factor(Temp_class)23  0.54244    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mass_mod_noMonth)      # residuals plot

ranova(mass_mod_noMonth)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + Month:as.factor(Temp_class)
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>          15 331.56 -633.11                         
## (1 | ID_fish)   14 115.72 -203.43 431.68  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mass model with random intercepts by temp treatment
mass_mod_trtint <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (0 + Month_int|Temp_class), 
                 data=ds1, REML=T)

summary(mass_mod_trtint)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (0 +  
##     Month_int | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -231.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.77281 -0.70174  0.00548  0.69951  2.82668 
## 
## Random effects:
##  Groups     Name      Variance Std.Dev.
##  Temp_class Month_int 0.12294  0.3506  
##  Residual             0.03292  0.1814  
## Number of obs: 509, groups:  Temp_class, 2
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                             2.277e-01  3.554e-01  4.960e+02   0.641
## MonthJune                               1.658e-01  3.520e-01  4.960e+02   0.471
## MonthJuly                               2.258e-01  7.020e-01  4.960e+02   0.322
## MonthAugust                             1.583e-01  1.052e+00  4.960e+02   0.150
## MonthSeptember                          3.704e-01  1.403e+00  4.960e+02   0.264
## as.factor(Temp_class)23                -7.388e-02  4.971e-01  4.960e+02  -0.149
## SexF                                    1.055e-01  2.758e-02  4.960e+02   3.824
## SexM                                    3.331e-02  2.858e-02  4.960e+02   1.165
## Density                                 5.294e-03  6.385e-03  4.960e+02   0.829
## MonthJune:as.factor(Temp_class)23      -2.945e-02  4.985e-01  4.960e+02  -0.059
## MonthJuly:as.factor(Temp_class)23      -4.927e-04  9.930e-01  4.960e+02   0.000
## MonthAugust:as.factor(Temp_class)23     7.486e-02  1.489e+00  4.960e+02   0.050
## MonthSeptember:as.factor(Temp_class)23 -1.122e-02  1.984e+00  4.960e+02  -0.006
##                                        Pr(>|t|)    
## (Intercept)                            0.521955    
## MonthJune                              0.637892    
## MonthJuly                              0.747808    
## MonthAugust                            0.880518    
## MonthSeptember                         0.791861    
## as.factor(Temp_class)23                0.881916    
## SexF                                   0.000148 ***
## SexM                                   0.244381    
## Density                                0.407426    
## MonthJune:as.factor(Temp_class)23      0.952906    
## MonthJuly:as.factor(Temp_class)23      0.999604    
## MonthAugust:as.factor(Temp_class)23    0.959909    
## MonthSeptember:as.factor(Temp_class)23 0.995491    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mass_mod_trtint)      # residuals plot

ranova(mass_mod_trtint)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (0 + Month_int | Temp_class) + Month:as.factor(Temp_class)
##                                           npar logLik     AIC        LRT Df
## <none>                                      15 115.72 -201.43              
## Month_int in (0 + Month_int | Temp_class)   15 115.72 -201.43 5.0591e-12  0
##                                           Pr(>Chisq)
## <none>                                              
## Month_int in (0 + Month_int | Temp_class)

Plot predicted mean level trend in mass over time, by temperature

pred <- ggpredict(mass_mod, terms = c("Month", "Temp_class")) # generate trend lines for the two temperature treatments

ggplot() +
    geom_errorbar(
      data = pred, 
      aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add 95% confidence interval bars
  geom_point(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group)) +                                    # predicted level by temp*month
  geom_line(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group)) +                                    # show level trend over time
  labs(
    x = "Month", 
    y = "Predicted log_mass", 
    color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00"))

F test on interaction effect

anova(mass_mod, type = "3") # type 3 anova because we expect a meaningful interaction
## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Month                       1.80879 0.45220     4 215.74 129.7554 < 2.2e-16 ***
## as.factor(Temp_class)       0.01454 0.01454     1 107.54   4.1733   0.04351 *  
## Sex                         0.03252 0.01626     2 110.81   4.6664   0.01133 *  
## Density                     0.01128 0.01128     1 468.64   3.2366   0.07265 .  
## Month:as.factor(Temp_class) 0.10964 0.02741     4 212.06   7.8649 6.274e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Generate fitted data

z_mass<-augment(mass_mod)                             # get variables and fitted values from model
head(z_mass)
## # A tibble: 6 × 19
##   .rownames log_mass Month `as.factor(Temp_class)` Sex   Density Month_int
##   <chr>        <dbl> <fct> <fct>                   <fct>   <dbl>     <dbl>
## 1 1           0.423  May   23                      F           8         1
## 2 2          -0.0915 May   23                      M           8         1
## 3 3           0.250  May   23                      F           8         1
## 4 4           0.436  May   23                      M           8         1
## 5 5           0.360  May   23                      M           8         1
## 6 7           0.310  May   23                      M           8         1
## # ℹ 12 more variables: ID_fish <dbl>, .fitted <dbl>, .resid <dbl>, .hat <dbl>,
## #   .cooksd <dbl>, .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>,
## #   .sqrtrwt <dbl>, .weights <dbl>, .wtres <dbl>
##### check assumptions: plot all residual values against the fitted 
ggplot(z_mass, aes(x = .fitted, y = .resid)) +   # same plot as before
  geom_point(size = 2) 

FIGURE 1

fig1.1 <- ggplot(z_mass, aes(x = Month, y = log_mass, group = ID_fish)) + 
  geom_line(aes(y = .fitted, color = `as.factor(Temp_class)`),             # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  labs(x = "Month", y = expression(log[10]~(mass~(g))), color = expression("Temperature ("*degree*C*")")) +
  ylim(-0.24,1) +
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

fig1.2 <- ggplot() +
 geom_errorbar(
      data = pred, 
      aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2, alpha = 0.3) +      # add 95% confidence interval bars
  geom_line(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  labs(
    x = "Month", 
    y = expression(log[10]~(mass~(g))), 
    color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  ylim(-0.24,1) +
  scale_color_manual(values = c("#78B7C5", "#F21A00"))

fig1 <- ggarrange(print(fig1.1 + rremove("xlab") + rremove ("ylab")), 
                  print(fig1.2 + rremove("xlab") + rremove ("ylab")), 
                           ncol = 2, nrow = 1, widths = c(1, 1.1),
                           align = "v", 
                           common.legend = TRUE, 
                           legend = "right")

require(grid)

jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
  fig1, 
  bottom = textGrob("Month"),
  left = textGrob(
    expression(log[10]~(mass~(g))), 
    rot = 90,          # Rotate the y-axis label 90 degrees
    gp = gpar(fontsize = 12)  # Optionally adjust the font size
  )
)

setwd(figures_wd)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
  fig1, 
  bottom = textGrob("Month"),
  left = textGrob(
    expression(log[10]~(mass~(g))), 
    rot = 90,          # Rotate the y-axis label 90 degrees
    gp = gpar(fontsize = 12)  # Optionally adjust the font size
  )
)
dev.off()
## png 
##   2

SMR MODEL

Visualize SMR data

hist(ds1$SMR)        # plot all raw SMR values in a histogram, note skewness

hist(log10(ds1$SMR)) # plot all log-SMR values in a histogram

ggplot(ds1, aes(x = mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against mass
  geom_point(size = 2) +                                        # size of data points
  facet_wrap(. ~ Temp_class) +                                  # produces one plot per temperature treatment
  theme_bw() 

Visualize log-mass and SMR relationship

ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against log10(mass)
  geom_point(size = 2) +                                            # size of data points
  facet_wrap(. ~ Temp_class) +                                      # produces one plot per temperature treatment
  theme_bw() 

Note that this relationship is non-linear, we expect this to be described by a power function y = a*x^b

This reflects our understanding of the relationship between mass and metabolism, which can be described by the power function y = a*x^b

ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against log10(mass)
  geom_point(size = 2) +                                            # size of data points
  facet_wrap(. ~ Temp_class) +                                      # produces one plot per temperature treatment
  theme_bw() +
  scale_y_log10()                                                   # plot data with log10 scaled y axis

The trends become linear when SMR is presented on log scale

Univariate SMR models for two temperature treatments

# model structure using log10 SMR and mass for fish at 18C only
SMR18_mod<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes by month, random intercepts by ID_fish
                 data=ds18, REML = T)

summary(SMR18_mod)                                               # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -553.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6860 -0.7401  0.0462  0.7148  2.7433 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  ID_fish  (Intercept) 1.063e-17 3.261e-09
##  Residual             9.832e-03 9.916e-02
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.67878    0.01415 315.00000  -47.96   <2e-16 ***
## log_mass      0.78058    0.02562 315.00000   30.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(SMR18_mod)                                                  # plot model residuals

ranova(SMR18_mod)                                                # anova-like test shows insignificant random slopes
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC LRT Df Pr(>Chisq)
## <none>           4 276.78 -545.56                  
## (1 | ID_fish)    3 276.78 -547.56   0  1          1
SMR18_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish),          # random slopes term removed
                 data=ds18, REML = T)

summary(SMR18_mod_1)                                             # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -553.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6860 -0.7401  0.0462  0.7148  2.7433 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  ID_fish  (Intercept) 1.063e-17 3.261e-09
##  Residual             9.832e-03 9.916e-02
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.67878    0.01415 315.00000  -47.96   <2e-16 ***
## log_mass      0.78058    0.02562 315.00000   30.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(SMR18_mod_1)                                                # plot model residuals

ranova(SMR18_mod_1)                                              # anova-like test shows insignificant random intercepts, but retain in model for testing
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC LRT Df Pr(>Chisq)
## <none>           4 276.78 -545.56                  
## (1 | ID_fish)    3 276.78 -547.56   0  1          1
# model structure using log10 SMR and mass for fish at 23C only
SMR23_mod<-lmer(log10(SMR) ~ log_mass + (1 + Month_int|ID_fish), # random slopes by month, random intercepts by ID_fish
                 data=ds23, REML=T)

summary(SMR23_mod)                                               # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -649.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1396 -0.5503 -0.0023  0.5686  3.2950 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 0.0052305 0.07232       
##           Month_int   0.0003171 0.01781  -0.94
##  Residual             0.0056041 0.07486       
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.47671    0.01133 101.03860  -42.07   <2e-16 ***
## log_mass      0.61389    0.02319 128.47109   26.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.876
plot(SMR23_mod)                                                  # plot model residuals

ranova(SMR23_mod)                                                # anova-like test shows insignificant random slopes
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
##                                        npar logLik     AIC    LRT Df Pr(>Chisq)
## <none>                                    6 324.70 -637.41                     
## Month_int in (1 + Month_int | ID_fish)    4 320.96 -633.92 7.4863  2    0.02368
##                                         
## <none>                                  
## Month_int in (1 + Month_int | ID_fish) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SMR23_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish),          # random slopes term removed
                 data=ds23, REML=T)

summary(SMR23_mod_1)                                             # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -641.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0810 -0.5308 -0.0445  0.5962  3.1721 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0006479 0.02545 
##  Residual             0.0064034 0.08002 
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.47203    0.01064 133.59553  -44.38   <2e-16 ***
## log_mass      0.60664    0.02244 178.85157   27.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.851
plot(SMR23_mod_1)                                                # plot model residuals

ranova(SMR23_mod_1)                                              # anova-like test shows insignificant random intercepts, but retain in model for testing
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC    LRT Df Pr(>Chisq)  
## <none>           4 320.96 -633.92                       
## (1 | ID_fish)    3 318.78 -631.56 4.3593  1    0.03681 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Combined SMR model

# model using log10 SMR and the interaction of log10 mass and temperature treatment
SMR_mod <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|ID_fish),  
                 data=ds1, REML=T)
plot(SMR_mod)                                             # generate model output

summary(SMR_mod)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -927.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.82335 -0.64895  0.01139  0.65120  3.06128 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0002204 0.01485 
##  Residual             0.0084623 0.09199 
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.648757   0.030786 191.298737 -21.073
## log_mass                           0.791230   0.025593 308.700979  30.916
## as.factor(Temp_class)23            0.178487   0.020589 235.734567   8.669
## SexF                              -0.025831   0.014616 181.956942  -1.767
## SexM                              -0.013551   0.014938 167.180010  -0.907
## Density                           -0.002235   0.003289 173.189264  -0.680
## log_mass:as.factor(Temp_class)23  -0.126302   0.039605 277.914923  -3.189
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          7.13e-16 ***
## SexF                              0.07886 .  
## SexM                              0.36562    
## Density                           0.49769    
## log_mass:as.factor(Temp_class)23  0.00159 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.449                                   
## as.fc(T_)23 -0.379  0.622                            
## SexF        -0.392 -0.194 -0.067                     
## SexM        -0.413 -0.066 -0.014  0.790              
## Density     -0.834  0.120  0.120  0.123  0.109       
## lg_:.(T_)23  0.340 -0.637 -0.898  0.046 -0.038 -0.105
ranova(SMR_mod)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
##               npar logLik     AIC     LRT Df Pr(>Chisq)
## <none>           9 463.94 -909.88                      
## (1 | ID_fish)    8 463.67 -911.34 0.54165  1     0.4618
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
SMR_mod_month <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(SMR_mod_month)                                             # generate model output

summary(SMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -929.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.85707 -0.62933  0.03354  0.64009  3.03111 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 1.442e-03 0.037971      
##           Month_int   4.593e-05 0.006777 -1.00
##  Residual             8.259e-03 0.090880      
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.644790   0.031143 232.514432 -20.704
## log_mass                           0.793677   0.025854 340.474164  30.699
## as.factor(Temp_class)23            0.178273   0.021153 199.523001   8.428
## SexF                              -0.026653   0.015115 151.602600  -1.763
## SexM                              -0.015453   0.015411 150.798427  -1.003
## Density                           -0.002854   0.003294 238.799373  -0.866
## log_mass:as.factor(Temp_class)23  -0.128338   0.039994 319.841634  -3.209
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          6.95e-15 ***
## SexF                              0.07985 .  
## SexM                              0.31760    
## Density                           0.38713    
## log_mass:as.factor(Temp_class)23  0.00147 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.452                                   
## as.fc(T_)23 -0.379  0.622                            
## SexF        -0.404 -0.197 -0.065                     
## SexM        -0.427 -0.069 -0.008  0.799              
## Density     -0.826  0.120  0.115  0.128  0.114       
## lg_:.(T_)23  0.345 -0.638 -0.902  0.047 -0.040 -0.107
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
ranova(SMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC    LRT Df Pr(>Chisq)
## <none>                                   11 464.62 -907.24                     
## Month_int in (1 + Month_int | ID_fish)    9 463.94 -909.88 1.3653  2     0.5053
# model with random slopes by temp treatment
SMR_mod_trt <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0+Temp_class|ID_fish),  
                 data=ds1, REML=T)

plot(SMR_mod_trt)                                             # generate model output

summary(SMR_mod_trt)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -927.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83458 -0.62750  0.01992  0.65194  3.04447 
## 
## Random effects:
##  Groups   Name         Variance  Std.Dev. Corr
##  ID_fish  Temp_class18 0.0001802 0.01342      
##           Temp_class23 0.0002780 0.01667  0.53
##  Residual              0.0084649 0.09201      
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.648557   0.030753 188.423517 -21.089
## log_mass                           0.791050   0.025482 219.289631  31.044
## as.factor(Temp_class)23            0.178657   0.020645 208.685088   8.654
## SexF                              -0.025952   0.014606 182.703141  -1.777
## SexM                              -0.013653   0.014928 167.741231  -0.915
## Density                           -0.002239   0.003292 171.722340  -0.680
## log_mass:as.factor(Temp_class)23  -0.126667   0.039728 233.094125  -3.188
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          1.35e-15 ***
## SexF                              0.07727 .  
## SexM                              0.36170    
## Density                           0.49730    
## log_mass:as.factor(Temp_class)23  0.00163 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.448                                   
## as.fc(T_)23 -0.374  0.617                            
## SexF        -0.391 -0.196 -0.067                     
## SexM        -0.413 -0.065 -0.013  0.789              
## Density     -0.835  0.120  0.119  0.122  0.108       
## lg_:.(T_)23  0.336 -0.632 -0.898  0.047 -0.037 -0.104
ranova(SMR_mod_trt)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC      LRT Df
## <none>                                     11 463.95 -905.91            
## Temp_class in (0 + Temp_class | ID_fish)    9 463.94 -909.88 0.028652  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)     0.9858
# model with random intercepts by temp treatment
SMR_mod_trtint <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)

plot(SMR_mod_trtint)                                             # generate model output

summary(SMR_mod_trtint)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -927.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86266 -0.65384  0.00029  0.65027  3.09859 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  Temp_class (Intercept) 0.0002948 0.01717 
##  Residual               0.0086736 0.09313 
## Number of obs: 509, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                      -6.473e-01  3.448e-02  1.288e-08 -18.773
## log_mass                          7.897e-01  2.513e-02  5.020e+02  31.428
## as.factor(Temp_class)23           1.768e-01  3.152e-02  2.247e-09   5.609
## SexF                             -2.541e-02  1.419e-02  5.020e+02  -1.791
## SexM                             -1.322e-02  1.447e-02  5.020e+02  -0.913
## Density                          -2.402e-03  3.186e-03  5.020e+02  -0.754
## log_mass:as.factor(Temp_class)23 -1.225e-01  3.879e-02  5.020e+02  -3.158
##                                  Pr(>|t|)    
## (Intercept)                       1.00000    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           1.00000    
## SexF                              0.07398 .  
## SexM                              0.36158    
## Density                           0.45117    
## log_mass:as.factor(Temp_class)23  0.00168 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.390                                   
## as.fc(T_)23 -0.483  0.399                            
## SexF        -0.343 -0.194 -0.041                     
## SexM        -0.363 -0.065 -0.006  0.797              
## Density     -0.720  0.114  0.076  0.126  0.111       
## lg_:.(T_)23  0.298 -0.639 -0.576  0.045 -0.039 -0.105
ranova(SMR_mod_trtint)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik     AIC        LRT Df Pr(>Chisq)
## <none>              9 463.67 -909.34                         
## (1 | Temp_class)    8 463.67 -911.34 1.1369e-13  1          1

Again, our ranova analysis shows that the random intercepts by ID_fish are not significant to this model. However, this will be retained because repeated measures were made on the same individuals, and so that the model can be included in the below covariance matrix.

Plot SMR vs. log mass mean predicted level

pred_SMR <- ggpredict(SMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass", y = "Predicted SMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

F test on interaction effect

anova(SMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                 Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)
## log_mass                       11.1822 11.1822     1 274.32 1321.4079 < 2.2e-16
## as.factor(Temp_class)           0.6359  0.6359     1 235.74   75.1494 7.134e-16
## Sex                             0.0318  0.0159     2 125.25    1.8793  0.156984
## Density                         0.0039  0.0039     1 173.19    0.4618  0.497689
## log_mass:as.factor(Temp_class)  0.0861  0.0861     1 277.92   10.1700  0.001591
##                                   
## log_mass                       ***
## as.factor(Temp_class)          ***
## Sex                               
## Density                           
## log_mass:as.factor(Temp_class) ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Generate fitted data

z_SMR<-augment(SMR_mod)                             # get variables and fitted values from model
head(z_SMR)
## # A tibble: 6 × 18
##   .rownames `log10(SMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1              -0.0733   0.423  23                      F           8      81
## 2 2              -0.546   -0.0915 23                      M           8      82
## 3 3              -0.276    0.250  23                      F           8      83
## 4 4              -0.0948   0.436  23                      M           8      84
## 5 5              -0.151    0.360  23                      M           8      85
## 6 7              -0.296    0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_SMR <- z_SMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )
##### check assumptions: plot all residual values against the fitted 
ggplot(z_SMR, aes(x = .fitted, y = .resid)) +   # same plot as before
  geom_point(size = 2) 

ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10()+
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

At this stage, we have found a strong model for SMR which captures variance by treatment and uses centered mass. The next step is to repeat this modeling approach for RMR and MMR.

RMR MODEL

# model using log10 RMR and the interaction of log10 mass and temperature treatment
RMR_mod <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod)

summary(RMR_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1011.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5876 -0.6152  0.0310  0.6508  2.7643 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.001429 0.03781 
##  Residual             0.006328 0.07955 
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.446375   0.032520 245.248665 -13.726
## log_mass                           0.714771   0.025398 410.246000  28.142
## as.factor(Temp_class)23            0.152471   0.021214 293.807345   7.187
## SexF                              -0.017279   0.015909 157.912118  -1.086
## SexM                              -0.016549   0.016440 144.900150  -1.007
## Density                           -0.003568   0.003490 252.616585  -1.022
## log_mass:as.factor(Temp_class)23  -0.138010   0.039679 393.028274  -3.478
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          5.48e-12 ***
## SexF                             0.279089    
## SexM                             0.315800    
## Density                          0.307526    
## log_mass:as.factor(Temp_class)23 0.000561 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.449                                   
## as.fc(T_)23 -0.360  0.595                            
## SexF        -0.385 -0.183 -0.073                     
## SexM        -0.396 -0.067 -0.036  0.756              
## Density     -0.843  0.155  0.123  0.106  0.095       
## lg_:.(T_)23  0.321 -0.630 -0.866  0.043 -0.033 -0.108
ranova(RMR_mod)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>           9 505.67 -993.33                         
## (1 | ID_fish)    8 494.08 -972.16 23.174  1   1.48e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
RMR_mod_month <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod_month)                                             # generate model output

summary(RMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1011.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6009 -0.6090  0.0215  0.6603  2.7544 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 1.876e-03 0.043312      
##           Month_int   3.614e-06 0.001901 -1.00
##  Residual             6.317e-03 0.079482      
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.445366   0.032513 245.812680 -13.698
## log_mass                           0.716091   0.025404 410.454083  28.188
## as.factor(Temp_class)23            0.150874   0.021468 219.587640   7.028
## SexF                              -0.017492   0.016095 137.340763  -1.087
## SexM                              -0.017322   0.016597 133.136674  -1.044
## Density                           -0.003705   0.003468 238.248975  -1.068
## log_mass:as.factor(Temp_class)23  -0.137589   0.039653 393.788850  -3.470
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          2.61e-11 ***
## SexF                             0.279008    
## SexM                             0.298544    
## Density                          0.286492    
## log_mass:as.factor(Temp_class)23 0.000578 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.452                                   
## as.fc(T_)23 -0.360  0.596                            
## SexF        -0.392 -0.187 -0.073                     
## SexM        -0.405 -0.070 -0.032  0.763              
## Density     -0.838  0.156  0.120  0.110  0.099       
## lg_:.(T_)23  0.325 -0.630 -0.870  0.045 -0.033 -0.110
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
ranova(RMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC     LRT Df
## <none>                                   11 505.81 -989.62           
## Month_int in (1 + Month_int | ID_fish)    9 505.67 -993.33 0.28935  2
##                                        Pr(>Chisq)
## <none>                                           
## Month_int in (1 + Month_int | ID_fish)     0.8653
# model with random intercepts by temp treatment
RMR_mod_trtint <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)
plot(RMR_mod_trtint)

summary(RMR_mod_trtint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -988.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9116 -0.6218  0.0056  0.6363  3.1738 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Temp_class (Intercept) 0.000744 0.02728 
##  Residual               0.007684 0.08766 
## Number of obs: 509, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                      -4.465e-01  3.919e-02  1.691e-08 -11.394
## log_mass                          7.095e-01  2.365e-02  5.020e+02  30.001
## as.factor(Temp_class)23           1.387e-01  4.296e-02  6.104e-09   3.229
## SexF                             -1.615e-02  1.336e-02  5.020e+02  -1.209
## SexM                             -1.567e-02  1.362e-02  5.020e+02  -1.150
## Density                          -3.326e-03  2.999e-03  5.020e+02  -1.109
## log_mass:as.factor(Temp_class)23 -1.078e-01  3.651e-02  5.020e+02  -2.954
##                                  Pr(>|t|)    
## (Intercept)                       1.00000    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           1.00000    
## SexF                              0.22718    
## SexM                              0.25062    
## Density                           0.26781    
## log_mass:as.factor(Temp_class)23  0.00328 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.323                                   
## as.fc(T_)23 -0.563  0.276                            
## SexF        -0.284 -0.194 -0.029                     
## SexM        -0.300 -0.065 -0.004  0.797              
## Density     -0.596  0.114  0.052  0.126  0.111       
## lg_:.(T_)23  0.247 -0.639 -0.397  0.045 -0.039 -0.105
ranova(RMR_mod_trtint)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik     AIC        LRT Df Pr(>Chisq)
## <none>              9 494.08 -970.16                         
## (1 | Temp_class)    8 494.08 -972.16 1.1369e-13  1          1
# model with random slopes by temp treatment
RMR_mod_trt <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod_trt)

summary(RMR_mod_trt)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1011.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5461 -0.6075  0.0311  0.6356  2.8170 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr
##  ID_fish  Temp_class18 0.001623 0.04028      
##           Temp_class23 0.001110 0.03331  0.06
##  Residual              0.006328 0.07955      
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.445596   0.032521 243.431641 -13.702
## log_mass                           0.715638   0.025664 339.122224  27.884
## as.factor(Temp_class)23            0.150977   0.020968 254.500115   7.200
## SexF                              -0.017989   0.015932 154.494731  -1.129
## SexM                              -0.017352   0.016446 142.933042  -1.055
## Density                           -0.003650   0.003464 233.113869  -1.054
## log_mass:as.factor(Temp_class)23  -0.134603   0.039235 317.874051  -3.431
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          6.76e-12 ***
## SexF                             0.260576    
## SexM                             0.293156    
## Density                          0.293117    
## log_mass:as.factor(Temp_class)23 0.000682 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.453                                   
## as.fc(T_)23 -0.375  0.609                            
## SexF        -0.390 -0.175 -0.072                     
## SexM        -0.398 -0.067 -0.037  0.762              
## Density     -0.839  0.154  0.128  0.111  0.098       
## lg_:.(T_)23  0.334 -0.644 -0.869  0.036 -0.035 -0.114
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
ranova(RMR_mod_trt)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC     LRT Df
## <none>                                     11 505.87 -989.74           
## Temp_class in (0 + Temp_class | ID_fish)    9 505.67 -993.33 0.41398  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)      0.813
# plot RMR vs. log mass mean predicted level 

pred_RMR <- ggpredict(RMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# plot the raw data as a cloud around the predictions

ggplot() +
  geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.3, position = position_jitter(width = 0.1)) +
  geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
             size = 3) +
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 2) +
  labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# F test on interaction effect

anova(RMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)
## log_mass                       6.5282  6.5282     1 386.99 1031.6186 < 2.2e-16
## as.factor(Temp_class)          0.3269  0.3269     1 293.81   51.6568 5.482e-12
## Sex                            0.0080  0.0040     2 120.77    0.6298 0.5344197
## Density                        0.0066  0.0066     1 252.62    1.0455 0.3075257
## log_mass:as.factor(Temp_class) 0.0766  0.0766     1 393.03   12.0973 0.0005614
##                                   
## log_mass                       ***
## as.factor(Temp_class)          ***
## Sex                               
## Density                           
## log_mass:as.factor(Temp_class) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# generate fitted data

z_RMR<-augment(RMR_mod)                                            # get variables and fitted values from model
head(z_RMR)
## # A tibble: 6 × 18
##   .rownames `log10(RMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1             -0.0212    0.423  23                      F           8      81
## 2 2             -0.386    -0.0915 23                      M           8      82
## 3 3             -0.147     0.250  23                      F           8      83
## 4 4              0.00674   0.436  23                      M           8      84
## 5 5             -0.0573    0.360  23                      M           8      85
## 6 7             -0.161     0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_RMR <- z_RMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )

# check assumptions: plot all residual values against the fitted 

ggplot(z_RMR, aes(x = .fitted, y = .resid)) +                      # same plot as before
  geom_point(size = 2) 

ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      
            alpha = 0.6) +
  theme_bw() +
  scale_y_log10() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

MMR MODEL

# model using log10 MMR and the interaction of log10 mass and temperature treatment

MMR_mod <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),  
                 data=ds1, REML=T)

plot(MMR_mod)

summary(MMR_mod) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1166.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5093 -0.5974  0.0398  0.5682  3.4530 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.001423 0.03772 
##  Residual             0.004466 0.06683 
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.013632   0.028910 251.274131  -0.472
## log_mass                           0.696766   0.022112 432.420918  31.511
## as.factor(Temp_class)23            0.047231   0.018743 292.436526   2.520
## SexF                              -0.015674   0.014420 139.112250  -1.087
## SexM                               0.008804   0.014953 127.598561   0.589
## Density                           -0.011717   0.003096 271.720060  -3.784
## log_mass:as.factor(Temp_class)23  -0.084810   0.034606 419.365154  -2.451
##                                  Pr(>|t|)    
## (Intercept)                       0.63766    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           0.01227 *  
## SexF                              0.27891    
## SexM                              0.55705    
## Density                           0.00019 ***
## log_mass:as.factor(Temp_class)23  0.01466 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.448                                   
## as.fc(T_)23 -0.354  0.585                            
## SexF        -0.387 -0.178 -0.074                     
## SexM        -0.395 -0.066 -0.042  0.747              
## Density     -0.843  0.166  0.123  0.101  0.090       
## lg_:.(T_)23  0.314 -0.628 -0.852  0.041 -0.032 -0.108
ranova(MMR_mod)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>           9 583.23 -1148.5                         
## (1 | ID_fish)    8 567.20 -1118.4 32.075  1  1.483e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
MMR_mod_month <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(MMR_mod_month)                                             # generate model output

summary(MMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1170.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1106 -0.5810  0.0524  0.5567  3.4994 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 0.0032626 0.05712       
##           Month_int   0.0001838 0.01356  -0.76
##  Residual             0.0040761 0.06384       
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.013991   0.029119 236.454925  -0.480
## log_mass                           0.694684   0.022919 329.612730  30.310
## as.factor(Temp_class)23            0.043012   0.019444 201.617723   2.212
## SexF                              -0.017277   0.014500 130.930308  -1.192
## SexM                               0.006459   0.014973 125.441099   0.431
## Density                           -0.011152   0.003128 218.545598  -3.566
## log_mass:as.factor(Temp_class)23  -0.077658   0.036235 275.968002  -2.143
##                                  Pr(>|t|)    
## (Intercept)                      0.631343    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          0.028082 *  
## SexF                             0.235604    
## SexM                             0.666953    
## Density                          0.000446 ***
## log_mass:as.factor(Temp_class)23 0.032973 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.433                                   
## as.fc(T_)23 -0.350  0.590                            
## SexF        -0.399 -0.182 -0.067                     
## SexM        -0.409 -0.065 -0.028  0.760              
## Density     -0.835  0.126  0.104  0.114  0.100       
## lg_:.(T_)23  0.311 -0.622 -0.870  0.038 -0.038 -0.091
ranova(MMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC  LRT Df Pr(>Chisq)
## <none>                                   11 585.47 -1148.9                   
## Month_int in (1 + Month_int | ID_fish)    9 583.23 -1148.5 4.47  2      0.107
# model with random intercepts by temp treatment
MMR_mod_trtint <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)

plot(MMR_mod_trtint)

summary(MMR_mod_trtint) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -1134.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2055 -0.6028  0.0431  0.6245  2.9068 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Temp_class (Intercept) 0.000000 0.00000 
##  Residual               0.005742 0.07578 
## Number of obs: 509, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.030099   0.024328 502.000000  -1.237
## log_mass                           0.699658   0.020445 502.000000  34.222
## as.factor(Temp_class)23            0.038596   0.016348 502.000000   2.361
## SexF                              -0.015362   0.011546 502.000000  -1.330
## SexM                               0.009194   0.011776 502.000000   0.781
## Density                           -0.009563   0.002592 502.000000  -3.689
## log_mass:as.factor(Temp_class)23  -0.065391   0.031558 502.000000  -2.072
##                                  Pr(>|t|)    
## (Intercept)                      0.216586    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          0.018612 *  
## SexF                             0.183961    
## SexM                             0.435357    
## Density                          0.000249 ***
## log_mass:as.factor(Temp_class)23 0.038765 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.449                                   
## as.fc(T_)23 -0.382  0.626                            
## SexF        -0.396 -0.194 -0.065                     
## SexM        -0.418 -0.065 -0.009  0.797              
## Density     -0.831  0.114  0.119  0.126  0.111       
## lg_:.(T_)23  0.344 -0.639 -0.903  0.045 -0.039 -0.105
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
ranova(MMR_mod_trtint)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik     AIC         LRT Df Pr(>Chisq)
## <none>              9  567.2 -1116.4                          
## (1 | Temp_class)    8  567.2 -1118.4 -2.2737e-13  1          1
# model with random slopes by temp treatment
MMR_mod_trt <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),  
                 data=ds1, REML=T)

plot(MMR_mod_trt)

summary(MMR_mod_trt) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1168.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3634 -0.6080  0.0343  0.5745  3.4158 
## 
## Random effects:
##  Groups   Name         Variance  Std.Dev. Corr 
##  ID_fish  Temp_class18 0.0018262 0.04273       
##           Temp_class23 0.0008205 0.02865  -0.12
##  Residual              0.0044550 0.06675       
## Number of obs: 509, groups:  ID_fish, 113
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.016069   0.028809 246.667733  -0.558
## log_mass                           0.698513   0.022560 378.875938  30.962
## as.factor(Temp_class)23            0.045074   0.018239 264.045593   2.471
## SexF                              -0.017409   0.014412 136.457536  -1.208
## SexM                               0.006586   0.014893 126.981343   0.442
## Density                           -0.011265   0.003029 231.897216  -3.719
## log_mass:as.factor(Temp_class)23  -0.079325   0.033707 333.250899  -2.353
##                                  Pr(>|t|)    
## (Intercept)                      0.577491    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          0.014095 *  
## SexF                             0.229152    
## SexM                             0.659105    
## Density                          0.000251 ***
## log_mass:as.factor(Temp_class)23 0.019184 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexF   SexM   Densty
## log_mass    -0.457                                   
## as.fc(T_)23 -0.390  0.615                            
## SexF        -0.400 -0.158 -0.068                     
## SexM        -0.401 -0.065 -0.044  0.762              
## Density     -0.833  0.162  0.134  0.111  0.098       
## lg_:.(T_)23  0.346 -0.660 -0.859  0.024 -0.039 -0.121
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
ranova(MMR_mod_trt)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC    LRT Df
## <none>                                     11 584.32 -1146.7          
## Temp_class in (0 + Temp_class | ID_fish)    9 583.23 -1148.5 2.1768  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)     0.3368
# plot MMR vs. log mass mean predicted level

pred_MMR <- ggpredict(MMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass", y = "Predicted MMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# Plot the raw data as a cloud around the predictions
ggplot() +
  geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.3, position = position_jitter(width = 0.1)) +
  geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
             size = 3) +
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 2) +
  labs(x = "log_mass_centre", y = "Predicted MMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# F test on interaction effect

anova(MMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)
## log_mass                       6.2104  6.2104     1 412.94 1390.6566 < 2.2e-16
## as.factor(Temp_class)          0.0284  0.0284     1 292.44    6.3497 0.0122725
## Sex                            0.0251  0.0126     2 108.54    2.8111 0.0645219
## Density                        0.0640  0.0640     1 271.72   14.3207 0.0001897
## log_mass:as.factor(Temp_class) 0.0268  0.0268     1 419.37    6.0062 0.0146630
##                                   
## log_mass                       ***
## as.factor(Temp_class)          *  
## Sex                            .  
## Density                        ***
## log_mass:as.factor(Temp_class) *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Generate fitted data

z_MMR<-augment(MMR_mod)                             # get variables and fitted values from model
head(z_MMR)
## # A tibble: 6 × 18
##   .rownames `log10(MMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1               0.154    0.423  23                      F           8      81
## 2 2               0.0458  -0.0915 23                      M           8      82
## 3 3               0.159    0.250  23                      F           8      83
## 4 4               0.215    0.436  23                      M           8      84
## 5 5               0.227    0.360  23                      M           8      85
## 6 7               0.0919   0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_MMR <- z_MMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )

# check assumptions: plot all residual values against the fitted 
ggplot(z_MMR, aes(x = .fitted, y = .resid)) +       # same plot as before
  geom_point(size = 2) 

ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),            
            alpha = 0.6) +
  theme_bw() +
  scale_y_log10()+
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

FIGURE 2

Generate a figure that displays mass-scaling relationships for all three metabolic rates at both temperatures

### FIG. 2 ###

fig2.1 <- ggplot() +
    geom_ribbon(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = SMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = "log10(mass (g))", y = expression(SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
 scale_y_log10(limits = c(0.15, 1.6))

fig2.2 <- ggplot() +
    geom_ribbon(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = "log10(mass (g))", y = expression(RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
  scale_y_log10(limits=c(0.25,2.1))

fig2.3 <- ggplot() +
    geom_ribbon(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
  scale_y_log10(limits=c(0.5, 5))

fig2.4 <- ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits = c(0.15, 1.6))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))


fig2.5 <- ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits=c(0.25,2.1))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

fig2.6 <- ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits=c(0.5, 5))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

library(ggpubr)

fig2 <- ggarrange(print(fig2.1 + rremove("xlab")), 
                  print(fig2.2 + rremove("xlab")), 
                  print(fig2.3 + rremove("xlab")),
                  print(fig2.4 + rremove("xlab") + rremove("ylab")),
                  print(fig2.5 + rremove("xlab") + rremove("ylab")),
                  print(fig2.6 + rremove("xlab") + rremove("ylab")),
                           ncol = 3, nrow = 2, widths = c(5, 5, 5),
                           align = "v", 
                           common.legend = TRUE, 
                           legend = "right")

require(grid)
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))

setwd(figures_wd)
jpeg(file = "Fig2.jpeg", width = 1440, height=960, units = "px")
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
dev.off()
## png 
##   2

UNIVARIATE MODEL OF FAT MASS

#subset ds1 to only include fish for which fat mass was measured
ds_fat <- ds1[which(!is.na(ds1$log_fatmass)),]

fat_mod <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + Sex + Density, data=ds_fat)

fat_mod1 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + SMR:Temp_class + Sex + Density, data=ds_fat)

fat_mod2 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + RMR:Temp_class + Sex + Density, data=ds_fat)

fat_mod3 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + MMR:Temp_class + Sex + Density, data=ds_fat)

summary(fat_mod)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83753 -0.09771  0.03729  0.12629  0.31222 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.714798   0.064714 -26.498  < 2e-16 ***
## Temp_class23                      -0.337119   0.062082  -5.430 9.31e-08 ***
## SexM                              -0.013164   0.018726  -0.703    0.482    
## Density                           -0.006225   0.006450  -0.965    0.335    
## Temp_class18:log_finalmass_centre  1.627573   0.065734  24.760  < 2e-16 ***
## Temp_class23:log_finalmass_centre  2.421200   0.083660  28.941  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1826 on 442 degrees of freedom
##   (115 observations deleted due to missingness)
## Multiple R-squared:  0.7806, Adjusted R-squared:  0.7781 
## F-statistic: 314.4 on 5 and 442 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

summary(fat_mod1)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     SMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82513 -0.08758  0.03480  0.12476  0.34453 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.709110   0.064990 -26.298  < 2e-16 ***
## Temp_class23                      -0.324979   0.062354  -5.212 2.88e-07 ***
## SexM                              -0.015560   0.018727  -0.831   0.4065    
## Density                           -0.006747   0.006456  -1.045   0.2966    
## Temp_class18:log_finalmass_centre  1.617051   0.075888  21.308  < 2e-16 ***
## Temp_class23:log_finalmass_centre  2.574283   0.114323  22.518  < 2e-16 ***
## Temp_class18:SMR                   0.008286   0.043880   0.189   0.8503    
## Temp_class23:SMR                  -0.136993   0.069860  -1.961   0.0505 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1822 on 440 degrees of freedom
##   (115 observations deleted due to missingness)
## Multiple R-squared:  0.7825, Adjusted R-squared:  0.779 
## F-statistic: 226.1 on 7 and 440 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

summary(fat_mod2)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     RMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82329 -0.08922  0.03455  0.12615  0.32535 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.7069333  0.0655146 -26.054  < 2e-16 ***
## Temp_class23                      -0.3273298  0.0627418  -5.217  2.8e-07 ***
## SexM                              -0.0162856  0.0188321  -0.865    0.388    
## Density                           -0.0068568  0.0064714  -1.060    0.290    
## Temp_class18:log_finalmass_centre  1.6226214  0.0766152  21.179  < 2e-16 ***
## Temp_class23:log_finalmass_centre  2.5482306  0.1180100  21.593  < 2e-16 ***
## Temp_class18:RMR                   0.0006086  0.0348560   0.017    0.986    
## Temp_class23:RMR                  -0.0865418  0.0567139  -1.526    0.128    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1825 on 440 degrees of freedom
##   (115 observations deleted due to missingness)
## Multiple R-squared:  0.7817, Adjusted R-squared:  0.7782 
## F-statistic: 225.1 on 7 and 440 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

summary(fat_mod3)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     MMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83208 -0.08588  0.03875  0.11949  0.32539 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.694523   0.066308 -25.555  < 2e-16 ***
## Temp_class23                      -0.338613   0.063490  -5.333 1.55e-07 ***
## SexM                              -0.013048   0.018730  -0.697    0.486    
## Density                           -0.007680   0.006523  -1.177    0.240    
## Temp_class18:log_finalmass_centre  1.679581   0.079663  21.084  < 2e-16 ***
## Temp_class23:log_finalmass_centre  2.489030   0.110451  22.535  < 2e-16 ***
## Temp_class18:MMR                  -0.020546   0.017446  -1.178    0.240    
## Temp_class23:MMR                  -0.024615   0.026652  -0.924    0.356    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1825 on 440 degrees of freedom
##   (115 observations deleted due to missingness)
## Multiple R-squared:  0.7817, Adjusted R-squared:  0.7782 
## F-statistic:   225 on 7 and 440 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

pred_fat <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Temp_class"))

# Plot the raw data
ggplot() +
  geom_point(data = ds_fat, aes(x = log_finalmass_centre, y = log_fatmass, group = factor(Temp_class), color = factor(Temp_class))) +
  geom_line(data = pred_fat, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) 

# F test on interaction term
anova(fat_mod)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                  Df Sum Sq Mean Sq  F value    Pr(>F)    
## Temp_class                        1  0.214  0.2144   6.4335   0.01154 *  
## Sex                               1  1.014  1.0144  30.4389 5.877e-08 ***
## Density                           1  0.009  0.0086   0.2594   0.61080    
## Temp_class:log_finalmass_centre   2 51.159 25.5793 767.5449 < 2.2e-16 ***
## Residuals                       442 14.730  0.0333                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod1)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                  Df Sum Sq Mean Sq  F value    Pr(>F)    
## Temp_class                        1  0.214  0.2144   6.4609   0.01137 *  
## Sex                               1  1.014  1.0144  30.5685 5.534e-08 ***
## Density                           1  0.009  0.0086   0.2605   0.61004    
## Temp_class:log_finalmass_centre   2 51.159 25.5793 770.8136 < 2.2e-16 ***
## Temp_class:SMR                    2  0.129  0.0644   1.9412   0.14476    
## Residuals                       440 14.601  0.0332                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod2)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                  Df Sum Sq Mean Sq  F value    Pr(>F)    
## Temp_class                        1  0.214  0.2144   6.4383   0.01151 *  
## Sex                               1  1.014  1.0144  30.4615 5.827e-08 ***
## Density                           1  0.009  0.0086   0.2596   0.61067    
## Temp_class:log_finalmass_centre   2 51.159 25.5793 768.1164 < 2.2e-16 ***
## Temp_class:RMR                    2  0.078  0.0388   1.1646   0.31302    
## Residuals                       440 14.653  0.0333                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod3)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                  Df Sum Sq Mean Sq  F value    Pr(>F)    
## Temp_class                        1  0.214  0.2144   6.4367   0.01152 *  
## Sex                               1  1.014  1.0144  30.4540 5.848e-08 ***
## Density                           1  0.009  0.0086   0.2595   0.61071    
## Temp_class:log_finalmass_centre   2 51.159 25.5793 767.9250 < 2.2e-16 ***
## Temp_class:MMR                    2  0.074  0.0370   1.1094   0.33066    
## Residuals                       440 14.656  0.0333                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This model is excluded from the following covariance matrix because we were only able to collect one measurement of fat mass from each fish, all at the very end of the experiment. Without an initial measurement of fat mass (logistically impossible as it was a lethal measurement), we are unable to assess within-subject residual correlation. For this reason, we have assessed it separately.

COVARIANCE MATRIX

build covariance matrix with brms

RN <- bf(mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred |a| ID_fish)) +
      bf(SMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
      bf(RMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
      bf(MMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
      set_rescor(T)

prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("lkj(2)", class = "cor"))

 fit1 <- brm(formula = RN, 
            data = ds1, 
            prior = prior1, 
            seed = 42, warmup = 500, iter = 10000, chains = 4, cores = 4,
            control=list(adapt_delta=0.97, max_treedepth = 15),
            save_ranef = T
           )

print(fit1, digits = 3)
##  Family: MV(gaussian, gaussian, gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred | a | ID_fish) 
##          SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##          RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##          MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##    Data: ds1 (Number of observations: 509) 
##   Draws: 4 chains, each with iter = 10000; warmup = 500; thin = 1;
##          total post-warmup draws = 38000
## 
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 113) 
##                                          Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept)                         0.180     0.013    0.156    0.207
## sd(mass1_Month_centred)                     0.042     0.004    0.035    0.051
## sd(SMR1_Intercept)                          0.015     0.007    0.002    0.028
## sd(RMR1_Intercept)                          0.033     0.005    0.024    0.043
## sd(MMR1_Intercept)                          0.037     0.005    0.027    0.047
## cor(mass1_Intercept,mass1_Month_centred)   -0.350     0.097   -0.529   -0.147
## cor(mass1_Intercept,SMR1_Intercept)        -0.300     0.273   -0.758    0.300
## cor(mass1_Month_centred,SMR1_Intercept)     0.434     0.243   -0.145    0.819
## cor(mass1_Intercept,RMR1_Intercept)        -0.055     0.154   -0.355    0.245
## cor(mass1_Month_centred,RMR1_Intercept)     0.252     0.142   -0.034    0.518
## cor(SMR1_Intercept,RMR1_Intercept)          0.512     0.256   -0.143    0.852
## cor(mass1_Intercept,MMR1_Intercept)        -0.090     0.150   -0.375    0.209
## cor(mass1_Month_centred,MMR1_Intercept)     0.311     0.135    0.034    0.561
## cor(SMR1_Intercept,MMR1_Intercept)          0.169     0.278   -0.416    0.669
## cor(RMR1_Intercept,MMR1_Intercept)          0.416     0.146    0.111    0.681
##                                           Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept)                      1.000     6949    13169
## sd(mass1_Month_centred)                  1.000     7814    16260
## sd(SMR1_Intercept)                       1.001     5973     8296
## sd(RMR1_Intercept)                       1.002     2822     6513
## sd(MMR1_Intercept)                       1.000    13058    18056
## cor(mass1_Intercept,mass1_Month_centred) 1.001     7166    14793
## cor(mass1_Intercept,SMR1_Intercept)      1.000    17582    20725
## cor(mass1_Month_centred,SMR1_Intercept)  1.000    20046    15852
## cor(mass1_Intercept,RMR1_Intercept)      1.000    14263    21929
## cor(mass1_Month_centred,RMR1_Intercept)  1.000    11969    21026
## cor(SMR1_Intercept,RMR1_Intercept)       1.003     1489     2408
## cor(mass1_Intercept,MMR1_Intercept)      1.000    18123    25293
## cor(mass1_Month_centred,MMR1_Intercept)  1.000    14145    21953
## cor(SMR1_Intercept,MMR1_Intercept)       1.001     1933     3914
## cor(RMR1_Intercept,MMR1_Intercept)       1.000     7547    17439
## 
## Regression Coefficients:
##                                          Estimate Est.Error l-95% CI u-95% CI
## mass1_Temp_class18                          0.175     0.061    0.056    0.294
## mass1_Temp_class23                          0.108     0.062   -0.014    0.231
## mass1_SexF                                  0.121     0.045    0.033    0.208
## mass1_SexM                                  0.050     0.048   -0.043    0.143
## mass1_Density                               0.010     0.006   -0.001    0.022
## mass1_Temp_class18:Month_centred_factor1    0.160     0.012    0.137    0.183
## mass1_Temp_class23:Month_centred_factor1    0.128     0.015    0.098    0.157
## mass1_Temp_class18:Month_centred_factor2    0.224     0.015    0.194    0.255
## mass1_Temp_class23:Month_centred_factor2    0.217     0.019    0.180    0.255
## mass1_Temp_class18:Month_centred_factor3    0.168     0.021    0.126    0.209
## mass1_Temp_class23:Month_centred_factor3    0.225     0.025    0.176    0.273
## mass1_Temp_class18:Month_centred_factor4    0.378     0.025    0.329    0.428
## mass1_Temp_class23:Month_centred_factor4    0.355     0.031    0.295    0.415
## SMR1_Temp_class18                          -0.845     0.036   -0.917   -0.775
## SMR1_Temp_class23                          -0.635     0.035   -0.704   -0.566
## SMR1_SexF                                  -0.028     0.015   -0.057    0.001
## SMR1_SexM                                  -0.015     0.015   -0.044    0.015
## SMR1_Density                               -0.002     0.003   -0.008    0.005
## SMR1_Temp_class18:log_mass_centre           0.803     0.031    0.743    0.865
## SMR1_Temp_class23:log_mass_centre           0.675     0.035    0.607    0.746
## RMR1_Temp_class18                          -0.607     0.035   -0.677   -0.538
## RMR1_Temp_class23                          -0.427     0.034   -0.495   -0.360
## RMR1_SexF                                  -0.017     0.015   -0.047    0.013
## RMR1_SexM                                  -0.016     0.016   -0.047    0.014
## RMR1_Density                               -0.005     0.003   -0.011    0.002
## RMR1_Temp_class18:log_mass_centre           0.714     0.029    0.657    0.770
## RMR1_Temp_class23:log_mass_centre           0.584     0.033    0.520    0.649
## MMR1_Temp_class18                          -0.188     0.033   -0.253   -0.122
## MMR1_Temp_class23                          -0.120     0.033   -0.185   -0.056
## MMR1_SexF                                  -0.017     0.014   -0.045    0.011
## MMR1_SexM                                   0.008     0.015   -0.021    0.037
## MMR1_Density                               -0.011     0.003   -0.017   -0.005
## MMR1_Temp_class18:log_mass_centre           0.705     0.027    0.652    0.758
## MMR1_Temp_class23:log_mass_centre           0.622     0.031    0.561    0.682
##                                           Rhat Bulk_ESS Tail_ESS
## mass1_Temp_class18                       1.000     5589    11325
## mass1_Temp_class23                       1.001     5162     9990
## mass1_SexF                               1.001     3826     7176
## mass1_SexM                               1.001     3950     7862
## mass1_Density                            1.000    11273    17979
## mass1_Temp_class18:Month_centred_factor1 1.000    12484    22473
## mass1_Temp_class23:Month_centred_factor1 1.000    13650    23584
## mass1_Temp_class18:Month_centred_factor2 1.001     7445    15542
## mass1_Temp_class23:Month_centred_factor2 1.000     7430    15925
## mass1_Temp_class18:Month_centred_factor3 1.001     6026    12963
## mass1_Temp_class23:Month_centred_factor3 1.001     6060    13153
## mass1_Temp_class18:Month_centred_factor4 1.001     5275    10983
## mass1_Temp_class23:Month_centred_factor4 1.001     5588    11585
## SMR1_Temp_class18                        1.000    11492    19639
## SMR1_Temp_class23                        1.000    11580    20326
## SMR1_SexF                                1.000    15446    23010
## SMR1_SexM                                1.000    16118    24116
## SMR1_Density                             1.000    18842    26708
## SMR1_Temp_class18:log_mass_centre        1.000    12958    21232
## SMR1_Temp_class23:log_mass_centre        1.000    14744    22700
## RMR1_Temp_class18                        1.000    10436    18121
## RMR1_Temp_class23                        1.001    10171    17456
## RMR1_SexF                                1.000    13066    21330
## RMR1_SexM                                1.000    13478    21505
## RMR1_Density                             1.000    16399    24852
## RMR1_Temp_class18:log_mass_centre        1.000    12740    20892
## RMR1_Temp_class23:log_mass_centre        1.000    14206    21255
## MMR1_Temp_class18                        1.000    14348    22135
## MMR1_Temp_class23                        1.000    14599    22742
## MMR1_SexF                                1.000    17061    25325
## MMR1_SexM                                1.000    17841    24465
## MMR1_Density                             1.000    20724    27685
## MMR1_Temp_class18:log_mass_centre        1.000    18451    25115
## MMR1_Temp_class23:log_mass_centre        1.000    21389    26997
## 
## Further Distributional Parameters:
##             Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma_mass1    0.060     0.003    0.055    0.065 1.000    18474    24645
## sigma_SMR1     0.092     0.003    0.086    0.098 1.000    16337    23784
## sigma_RMR1     0.081     0.003    0.075    0.086 1.000    11070    22175
## sigma_MMR1     0.068     0.002    0.063    0.073 1.000    24222    25735
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1)   -0.038     0.072   -0.180    0.104 1.000    15109    22918
## rescor(mass1,RMR1)   -0.066     0.069   -0.202    0.070 1.000    14908    23218
## rescor(SMR1,RMR1)     0.862     0.013    0.836    0.886 1.000    16327    25735
## rescor(mass1,MMR1)   -0.103     0.066   -0.231    0.029 1.000    18390    24837
## rescor(SMR1,MMR1)     0.169     0.047    0.075    0.259 1.000    21566    27134
## rescor(RMR1,MMR1)     0.238     0.047    0.145    0.327 1.000    22986    27361
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit1)

conditional_effects(fit1)

REPEATABILITY ESTIMATES

Use posterior distribution to generate repeatability estimates

View(posterior::summarise_draws(as_draws_array(fit1))) # view all predicted values of random effects

# random effects are fitted as standard dev, not variance, so we need to square them first
# var names to insert below: SMR, RMR, MMR, growth

# SMR repeatability

Vint_SMR <- (as_draws_array(fit1, variable = "sd_ID_fish__SMR1_Intercept"))^2
Vresid_SMR <- (as_draws_array(fit1, variable = "sigma_SMR1"))^2
R_SMR <- Vint_SMR / (Vint_SMR + Vresid_SMR)

"SMR"; mean(R_SMR); quantile(R_SMR, probs = c(0.025, 0.975))
## [1] "SMR"
## [1] 0.03092714
##         2.5%        97.5% 
## 0.0004320659 0.0890740509
# RMR repeatability

Vint_RMR <- (as_draws_array(fit1, variable = "sd_ID_fish__RMR1_Intercept"))^2
Vresid_RMR <- (as_draws_array(fit1, variable = "sigma_RMR1"))^2
R_RMR <- Vint_RMR / (Vint_RMR + Vresid_RMR)

"RMR"; mean(R_RMR); quantile(R_RMR, probs = c(0.025, 0.975))
## [1] "RMR"
## [1] 0.1458303
##       2.5%      97.5% 
## 0.07860133 0.23501854
# MMR repeatability

Vint_MMR <- (as_draws_array(fit1, variable = "sd_ID_fish__MMR1_Intercept"))^2
Vresid_MMR <- (as_draws_array(fit1, variable = "sigma_MMR1"))^2
R_MMR <- Vint_MMR / (Vint_MMR + Vresid_MMR)

"MMR"; mean(R_MMR); quantile(R_MMR, probs = c(0.025, 0.975))
## [1] "MMR"
## [1] 0.2306318
##      2.5%     97.5% 
## 0.1311474 0.3357818
# Growth rate repeatability
Vint_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Month_centred"))^2
Vslope_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Intercept"))^2
corr_growth <- (as_draws_array(fit1, variable = "cor_ID_fish__mass1_Intercept__mass1_Month_centred"))
Vresid_growth <- (as_draws_array(fit1, variable = "sigma_mass1"))^2
COVis_growth <- corr_growth*sqrt(Vslope_growth)*sqrt(Vint_growth)

x<- 0 #(can be 0 - 4)
VAR0 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R0_growth <- VAR0/(VAR0 + Vresid_growth)

"Growth - May"; mean(R0_growth); quantile(R0_growth, probs = c(0.025, 0.975))
## [1] "Growth - May"
## [1] 0.3336915
##      2.5%     97.5% 
## 0.2450998 0.4309421
x<- 1 #(can be 0 - 4)
VAR1 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R1_growth <- VAR1/(VAR1 + Vresid_growth)

"Growth - June"; mean(R1_growth); quantile(R1_growth, probs = c(0.025, 0.975))
## [1] "Growth - June"
## [1] 0.8876298
##      2.5%     97.5% 
## 0.8519065 0.9176981
x<- 2 #(can be 0 - 4)
VAR2 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R2_growth <- VAR2/(VAR2 + Vresid_growth)

"Growth - July"; mean(R2_growth); quantile(R2_growth, probs = c(0.025, 0.975))
## [1] "Growth - July"
## [1] 0.9705673
##      2.5%     97.5% 
## 0.9599677 0.9790658
x<- 3 #(can be 0 - 4)
VAR3 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R3_growth <- VAR3/(VAR3 + Vresid_growth)

"Growth - August"; mean(R3_growth); quantile(R3_growth, probs = c(0.025, 0.975))
## [1] "Growth - August"
## [1] 0.9869686
##      2.5%     97.5% 
## 0.9821577 0.9907833
x<- 4 #(can be 0 - 4)
VAR4 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R4_growth <- VAR4/(VAR4 + Vresid_growth)

"Growth - September"; mean(R4_growth); quantile(R4_growth, probs = c(0.025, 0.975))
## [1] "Growth - September"
## [1] 0.9927106
##      2.5%     97.5% 
## 0.9899955 0.9948520